Inhomogeneous atomic Bose-Fermi mixtures in cubic lattices 
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We determine the ground state properties of inhomogeneous mixtures of bosons and fermions in cubic lat- 
tices by studying the Bose-Fermi Hubbard model including parabolic confining potentials. We present the exact 
solution in the limit of vanishing hopping (ultradeep lattices) and study the resulting domain structure of com- 
posite particles. For finite hopping we determine the domain boundaries between Mott-insulator plateaux and 
hopping-dominated regions for lattices of arbitrary dimensionality within perturbation theory. The results are 
compared with a new numerical method that is based on a Gutzwiller variational approach for the bosons and an 
exact treatment for the fermions. The findings can be applied as a guideline for future experiments with trapped 
atomic Bose-Fermi mixtures in optical lattices. 
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Mixtures of ultracold bosonic and fermionic particles have 
attracted a considerable amount of attention in recent years, 
to a high extent triggered by the perspective of achieving 
prima facie transitions to superfluidity in systems of neutral 
fermionic atoms Spectacular progress has already been 
achieved in the experimental manipulation of cold atoms in 
optical lattices with an astonishingly high degree of control. 
Most prominently, the superfluid-Mott insulator phase transi- 
tion in systems of bosonic atoms in optical lattices has been 
experimentally observed Q], and the production of degenerate 
Fermi gases in optical lattices has been recently achieved 01 ■ 
A perspective of key interest in this field lies in the possibil- 
ity of discovering and probing new quantum phases of matter 
by combining ideas from the study of Bose-Fermi mixtures 
in solid state systems and of cold quantum gases in optical 
lattices USE* 00]. In Ref - 01' the Bose-Fermi Hubbard 
(BFH) Hamihonian has been introduced and derived from 
the underlying microscopic many-body Hamiltonian, linking 
the experimentally accessible quantities to the model parame- 
ters. A mean field argument has been presented for the onset 
of a bosonic superfluid transition, and in a numerical analy- 
sis the behavior of on-site quantities in several situations has 
been studied. In Refs. j5l 0] the phase diagram of homo- 
geneous boson-fermion mixtures in optical lattices has been 
studied in a mean field approach, and the existence of a com- 
plex structure of phases of composite fermionic particles has 
been suggested. In Ref. |7] stable supersolid phases have 
been predicted for homogeneous Bose-Fermi mixtures in two- 
dimensional lattices. Finally, Ref. yg] addresses the task of 
assessing the phase diagram of the BFH model using an exact 
diagonalization approach for systems of small size. 



The investigations in Refs. |5j |6l |j, |8|] are confined to the 
homogeneous case, i.e., to a translationally invariance BFH 
Hamiltonian. While this is a very reasonable approach as far 
as the discussion of the actual phases in the thermodynamical 
limit is concerned, it is not quite the one encountered experi- 
mentally, e.g. in the case of trapped, ultracold atomic gases. 
The external confining potential, superimposed to the optical 
lattice potential, breaks the translational symmetry, and leads 



to profound modifications of the phase structure by allowing 
for the appearance of spatial domains of coexisting different 
hases along the lattice, as rece ntly studied for pure bosonic 
and pure fermionic systems 1 10]. Studies of such inhomo- 
geneous systems are of immediate relevance for the interpre- 
tation of experimental findings, where some confinement in a 
trap is necessary. 

In the present paper we study the effects of an inhomo- 
geneous confining potential on Mott and superfluid regions 
emerging in systems of Bose-Fermi mixtures in regular lat- 
tices at zero temperature. We show that the model is exactly 
solvable in the limit of very strong lattices (vanishing bosonic 
and fermionic hopping), and analyze the related structure of 
domains of composite particles. We then consider the gen- 
eral case of finite hopping in Z?-dimensional lattices, study 
the bulk properties of the system in Landau theory and lo- 
cal density approximation (LDA), and determine the general 
phase boundaries of the different domains. Notably, we in- 
troduce a method to treat the bosons within a Gutzwiller-type 
ansatz fill fl2il and the fermions exactly, a versatile method 
that is applicable to several systems of this kind. This method 
allows us to present for the first time the domain structure of 
inhomogeneous atomic mixtures in confining potentials and 
the respective phase diagrams for the homogeneous case. 

Starting point of our analysis is the single-band BFH 
Hamiltonian foil , which captures the essential properties of 
dilute mixtures in optical lattices under fairly general as- 
sumptions on the tunable physical parameters |4j. The grand 
canonical BFH Hamiltonian reads 



H = 



-jb £ (s& + b]k) -j F j2 (/;./• + f}fi 

U BB »b(»b - !) + u bf £ 



B F 
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Here, b\ and fi are the on-site bosonic and fermionic an- 
nihilation operators, respectively, whereas fi l B = SJS, and 
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■hp = fjfi- Sites are associated with a cubic D-dimensional 
lattice with fixed spacing, and i = (ii,...,ip,) denotes a D- 
tuple labeling the coordinates of a site i with coordination 
number d (i.e., the number of nearest neighbors). The symbol 
denotes summation over pairs of nearest neighbors. The 
first two terms in Eq. Q describe independent bosonic and 
fermionic nearest-neighbor hopping with positive amplitudes 
Jb and Jp. The subsequent line represents on-site boson- 
boson and boson-fermion interactions. Finally, the first two 
terms of the last line incorporate the external confining po- 
tential, which, in typical experimental situations, can be taken 
to be harmonic. The origin of the lattice is chosen to be at 
the minimum of the trapping potential, assumed to be equal 
for bosons and fermions, so that V B = V F = Vi := Vo K | 2 - 
This Hamiltonian is a generalization to systems of bosons and 
fermions of the frequently employed Bose-Hubbard model ex- 
hibiting the Mott to superfluid phase transition in bosonic sys- 
tems fl2l f]~4L [I5I1 . Expressions linking the model parameters 
to quantities that can be tuned in an actual experimental sit- 
uation, such as the depth of the optical lattice and the atomic 
scattering lengths, are provided in Ref. |4j]. 

I) Exact solution with vanishing hopping. - A surprisingly 
rich situation is already encountered in the case of vanishing 
hopping: J B = Jf — 0. In this case the Hamiltonian Hq is 
simply a sum of single-site contributions, and the eigenstates 
of the BFH model are tensor products of number states with 
state vectors \ip) = \no, nx, ■ ■ ■ )\mo, mi, ■ ■ • )> where n, = 
0, 1, 2, ... and m, = 0, 1 represent the occupation number of 
bosons and fermions at site i, respectively. For simplicity of 
notation, we will fix the energy scale by setting U bb = 1- We 
have (ijj\H \ij}) = J2i( n i ~ n i + U B Fninii + V^n, + mi) - 
P-B^i — PF^i) ='■ E(ni,m,i), where for the ground state 
with state vector |-0o) the occupation numbers take the specific 
values 

= f max(0, [{l + p B -Vi)/2\), if E(fk,0) < E(n h l), 
1 \ max(0, [(1 + Pb -Vi- U B f)/2}), otherwise, 

= f 0, ifE(ni,0) < E(fn,l), 
1 I 1, otherwise, 

where [.] denotes the closest integer to the value in brackets. 
According to the above determination, several types of com- 
posite particles can be formed. Composites consisting of rhi 
fermions and hi bosons are formed at site i, see Fig.^ Con- 
nected domains with fixed integer particle number are formed 
and, depending on the interaction strength Ubf and the re- 
lation of the respective chemical potentials ps and pp, the 
fermions distribute around the center of the trap or are pushed 
outwards. 

II) Finite hopping: perturbative treatment. - We now turn 
to the strong coupling limit where also small but finite hop- 
ping is allowed for. In a wide range of physical parameters, 
the strength of the hopping for bosons and fermions are ap- 
proximately of the same value, for instance for neutral atoms 
in optical lattices 01 . We set subsequently J p = J 'b =■ J \ 
and treat the small positive parameter J as a perturbation. As 
in Refs. HI, we introduce a mean field approximation, which 
amounts to a replacement of the bosonic operator products 




FIG. 1 : Distribution of integer boson and fermion numbers for the 
case Jb — Jf = and Vo = 0.002 for a D-dimensional cubic 
lattice. This is encoded in the color as shown in the bar on the right 
hand side (number of bosons, number of fermions) as a function of 
the component i 1 , the chemical potential /j,b, and Ubf- For the left 
(right) figure, /j,f = 6/j.b/5 (hf = Ms/5) is chosen. 



in Eq. according to b\bj 1 — ► ip B *bj + b\il> B — ^> b ^b*> 
the complex numbers ip B being variational parameters mod- 
eling the influence of neighboring atoms with the physical in- 
terpretation of a superfluid parameter. We consider the re- 
sulting corrections to the ground state energy, (V'ol-^olV'o) = 
Y^i E(h i ,m i ), to second order in J. Moreover, to study 
bulk properties we will make use of the local density ap- 
proximation (LDA). This means taking for each lattice site 
ip B to be equal to the corresponding values at neighboring 
sites. This is well justified for a sufficiently shallow trap- 
ping potential. In this approximation, the ground state energy 
reads E = (ip \H \ipo) + &E B + AE F + 0{J 3 ), where 
AE B = 2JdJ2 i (|V>b| 2 (1 + Jdn)), with n = (4n, + 2c, + 
2)/( c ? — 1), Q = 1 — 2hi — Vi + pb — Ubf^, and d is the 
coordination number of a Z?-dimensional cubic lattice (d — 6 
in three dimensions). We are now in the position to apply the 
Landau argument to determine the phase boundaries within 
LDA. If 1 + Jdr B > 0, then the approximate energy func- 
tional is minimized by having \ip l B \ = 0, which corresponds 
to the incompressible Mott situation for the bosons. In turn, 
for 1 + Jdr l B < the minimization requires \t/j B | > 0, and the 
bosons are superfluid. Exploiting this property, we can deter- 
mine the phase boundary between the hopping-dominated and 
the Mott regime at each site, corresponding to J = — l/(dri). 
To find the boundaries for the fermions, we consider that for 
small J and within LDA the bosons alter the fermionic chem- 
ical potential, introducing an effective site dependent chem- 
ical potential a F = pp — Ubf^i — Vi. At each site i we 
then consider the corresponding (infinite) homogeneous prob- 
lem H' F = -JJ2(i,j)(fifi + fjfi)~ MfE;^f. which is 
appropriate for sufficiently shallow external potentials. This 
Hamiltonian is diagonal in Fourier space, so that the exact 



3 



0.1 
'd.05 


|0.M 
0.02 


<3M 
0.02 


|0.oi 

0.02 



Bosons 




Fermions 












1 








■ ■■■■ ' f 

1 


— ' -- c 








1 

;■ z 

■ illB- 


-«0 -50 -40 


-30 -20 | -10 


10 20 30 


40 



FIG. 2: The boundaries between Mott and hopping-dominated re- 
gions for D = 3 (d — 6) for both bosons (left) and fermions (right), 
as a function of the site index i\, and of the hopping J = Jb = Jf- 
The values of Ubf = 0.3 and ^f — (j.b/5 are chosen (corre- 
sponding to the topmost plot on the right of Fig. 0, and, from top 
to bottom, fi B = (0.7, 2, 2.4, 4, 6.8). In each plot, the J = axis 
corresponds to the plots of Fig.Q The white solid line depicts the 
phase boundaries as determined in section II. The dashed line re- 
produces the same plots, but for Ubf = 0. In the same diagram, 
the background color encodes the variance of the on-site density 
(j b /f = (((^s/f) 2 ) ~~ ("s/f) 2 ) 1 ^ 2 from the numerical variational 
analysis discussed in section III. Dark blue (gray) corresponds to the 
Mott region with a B / F = 0. 



spectrum is given by Sk = —pip — ^Jj2s=i cos(fca), where 
k = (ki, kjj), the lattice spacing being set to 1 without 
loss of generality. Therefore, when —p F — 2dJ > 0, the 
ground state has no fermions present, being obviously a Mott 
state. Similarly, for —p F + 2d J < the ground state is a Mott 
state with exactly one fermion at each site. Fig. [2] shows the 
phase regions for an inhomogeneous Bose-Fermi mixture in 
a three-dimensional lattice with a weakly confining parabolic 
potential. The solid lines depict the boundaries between Mott 
and hopping dominated regions, respectively, as evaluated us- 
ing the above approach. Not surprisingly, one observes that 
at the center of the trap, where the potential acquires its mini- 
mum, lower values of the hopping are needed for the transition 
to the hopping-dominated regime. For appropriate fixed J, 
different spatial domains develop from the center of the trap. 
Depending on the value of fig, one observes an alternating 
sequence of Mott and hopping-dominated domains. An im- 
portant new feature that emerges in inhomogeneous BFH sys- 
tems differing from the situation encountered in pure bosonic 
or fermionic systems is a modulation of the phase regions due 
to the boson-fermion interaction. This can be understood by 
comparing the phase boundaries for the interacting mixture 
with the non-interacting case Ubf = 0. The boundaries are 
represented as dashed lines in Fig. [2] For the chosen parame- 
ters, the presence of the fermions in the center of the trap is re- 
flected by a tendency to form Mott domains for bosons. Com- 
paring this functional behavior with the fermion number per 



site in the case of vanishing hopping as depicted in Fig.^we 
see that the state diagram for the bosons is modified when the 
fermion number per site is exactly one. In turn, the presence 
of the bosons heavily modifies the boundaries between the 
Mott and the hopping-dominated domains for the fermions: 
the hopping-dominated regions are pushed outwards, and the 
value of the integer boson occupation number per site in the 
Mott phase sets the scale of this phenomenon. 

Ill) Finite hopping: variational theory. - In Fig. [2] we have 
also represented the variance of the on-site densities cr B , F . 
They are determined using the following variational approach. 
We consider at each bulk site i the corresponding infinite ho- 
mogeneous lattice Hamiltonian, Hi. The minimization of 
{4>i\Hi\4>i) over all state vectors will be replaced by a min- 
imization over state vectors respecting the univalence super- 
selection rule, \4>i) = \4> B )\4>p)- For the bosonic sector we 
introduce a Gutzwiller-type ansatz, \<j) B ) = Yii S ni ^ni\ n i) 
(see, e.g., Refs. fill frill and references therein), where the 
b % n form a probability distribution at each site I, n; = 0, 1, ... . 
After an exact discrete Fourier transformation of the fermionic 
operators, // = -i afce lfcn , we have for each site i, 

(tilHilfr) = E B + ^4<4|4%|4), 

k 

(1,3) 

i i 

withe], = -4J J2s = i c °s( k s) - - U B F((j)B\n B \(l)B) - Vi. 
Therefore, the state vector \<j)%) = J\ k E « <0 ajjO) minimizes 
the energy expectation value at fixed Gutzwiller amplitudes, 

EUbi,b\r--) = E B + £ 4. 

To determine the ground state, we have to minimize E^ in at 
each site i. Because this energy functional is not convex, 
the energy landscape exhibits local minima and determining 
the ground state leads to a non-convex optimization problem. 
However, the problem can be solved numerically using a sim- 
ulated annealing method (73]. The regions with exactly van- 
ishing local variance, cr B , F , identify the respective Mott re- 
gions (dark blue/grey in Fig. [5}. Qualitatively, we obtain very 
similar results in the perturbative and in the variational treat- 
ments. The perturbative findings are valid for small hopping 
only, while the numerical analysis relies on the Gutzwiller 
ansatz for bosons, which is appropriate in high spatial dimen- 
sions (D = 3) and in the superfluid regime flU fl9il . The 
behavior shown in Fig.|2]is thus a genuine effect of the boson- 
fermion interaction in the mixture, as it is predicted by the 
considered approximation schemes. For a system of harmon- 
ically trapped bosons it has been shown that the appearance 
of a Mott-insulator domain within a shell of superfluid atoms 
leads to satellite peaks in the global momentum distribution 
l20ll . This feature is accessible in experiments and can in par- 
ticular be used as an indication for the effect of the fermions 
on the boson Mott transition. 
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FIG. 3: State diagram for central sites of bosons (left) and fermions 
(right) for Ubf = 0.3 and fiF = /Ub/5. Depicted are the phase 
boundaries as determined using the argument from section II in LDA 
(solid lines) and in the Gutzwiller variational theory (color encod- 
ing). The dashed lines correspond to boundaries for Ubf = 0, de- 
termined using the Landau argument from section II for the bosons, 
and exactly for the fermions. 

TV) Behavior at the center of the trap: bulk properties. - 
For the central sites, within LDA, the inhomogeneous case is 
equivalent to the homogeneous case. To interpret the find- 
ings, we first recall how the phase diagram for the fermions 
would look like in the homogeneous case in the absence of 
bosons. In this case the BFH model describes a spinless 
fermion system with hopping contributions only. It can be 
solved without approximation as before with the help of a dis- 
crete Fourier transformation. The Mott states with exactly 
one or zero fermions per site can be distinguished from the 
hopping-dominated states, yielding a linear behavior of the 
phase boundary as a function of J = Jp (see section II). This 
is depicted in Fig. [5] with a dashed line. Within the perturba- 
tive treatment, the effect of the bosons is to give rise to an ef- 
fective fermionic chemical potential, reflecting the change of 
the number of bosons per site in the Mott phase. This in turn 
leads to integer discontinuous jumps in the phase boundaries. 
In this way, the presence of the bosons modifies the fermionic 
phase diagram. In turn, the presence of the fermions modu- 
lates the phase diagram for the bosons as compared to the stan- 



dard mean-field phase diagram of the Bose-Hubbard-model. 
Notably, the lobes associated to different boson numbers per 
site in the Mott insulator do not necessarily touch the straight 
line corresponding to J = 0. Again, we have compared these 
findings with the results obtained from the numerical analysis 
introduced in section III. The general behavior of the regions 
with exactly vanishing density variance within Gutzwiller ap- 
proximation, and Mott regions according to the perturbative 
results is very similar. However, the discontinuities are less 
pronounced within the variational approximation. This is due 
to the fact that in perturbation theory the zeroth order contri- 
bution is manifestly discontinuous. We have compared this 
behavior with the results obtained from an exact diagonal- 
ization of the Hamiltonian for small systems l2lll . obtaining 
qualitatively identical conclusions. 

In conclusion, we have studied in detail the phase struc- 
ture of the ground state of trapped inhomogeneous Bose- 
Fermi mixtures in optical lattices. The inhomogeneity leads 
to domains of Mott plateaux and hopping-dominated regions, 
where a complex interplay between interacting bosons and 
fermions is displayed. Introducing a new numerical method 
that treats fermions without approximation, we were able to 
visualize for the first time the effects of this complex interplay 
on the domain structure of both species and present the phase 
diagrams for the homogeneous case. These results will be 
compared with DMRG-methods for one-dimensional lattices 
in forthcoming work. The findings reported in the present 
work should provide a guideline and should be amenable to 
direct testing in the upcoming experiments 01 with trapped 
mixtures of bosonic and fermionic atoms in optical lattices. 
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